function [Vnew,Tnew,Enew,ETnew,TEnew] = refine_mesh_local(V,T,E,ET,TE,edge_mark)
% function [Vnew,Tnew,Enew,ETnew,TEnew] = refine_mesh_local(V,T,E,ET,TE,edge_mark)
%
% adapt the mesh by refine edge indicated in flag vector edge_mark
%

% number of edges need to be refine
ne_mark = sum(edge_mark); 
if(ne_mark == 0) % no need todo refine
    Vnew = V;
    Tnew = T;
    Enew = E;
    ETnew = ET;
    TEnew = TE;
    return;   
end

n_eg_marked = edge_mark(TE(:,1)) + edge_mark(TE(:,2)) + edge_mark(TE(:,3));
n1 = sum(n_eg_marked == 1);
n2 = sum(n_eg_marked == 2);
n3 = sum(n_eg_marked == 3);

%% 1. we only treated with 1->2 and 1->4 cases, ignore n0 and n2 should be 0
% reserve memory spaces and initialize the mesh data
nT = size(T,1);  nV = size(V,1); nE = size(E,1);
Vnew =   [V;   zeros(ne_mark, 2)                               ];
Enew =   [E;    zeros(ne_mark + n1 + 2*n2 + 3*n3, 2)];
ETnew = [ET; zeros(ne_mark + n1 + 2*n2 + 3*n3, 2)];
Tnew =   [T;   zeros(n1 + 2*n2 + 3*n3, 3)                 ];
TEnew = [TE; zeros(n1 + 2*n2 + 3*n3, 3)                 ];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 2. generate new vertex on each flaged edges
eg_refine = find(edge_mark == 1);
gdxnew = ((nV + 1):(nV + ne_mark))';  % the global index for the new vertex
Vnew(gdxnew,:) = (V(E(eg_refine,1),:) + V(E(eg_refine,2),:))*0.5;% (V(E(eg_refine,1),:) + V(E(eg_refne,2),:))*0.5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 3. generate new edges
%
%                     v1
%                   /     \
%              e3/        \ e2
%                /           \
%            v2 --- e1 --- v3
% 
edge_cache = zeros(nE, 3);   % subedge infomation cache
% 1. add subedge 1 at the original place
Enew(eg_refine,:) = [E(eg_refine, 1) gdxnew];  
ETnew(eg_refine, 1) = 0;  ETnew(eg_refine,2) = 0;  % clear ETnew old info
edge_cache(eg_refine,1) = eg_refine; % replace edge with subedge1
% 2. add subedge 2 at the end
eg2_idx = ((nE+1):(nE+ne_mark))';
Enew(eg2_idx, :) = [E(eg_refine,2) gdxnew]; % fill edge 2 to list
edge_cache(eg_refine,2) = eg2_idx;       % cache the subedge 2 information

% index the center nodes
edge_cache(eg_refine,3) = gdxnew;      

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 4. generate new triangles, there are two cases:
if(n3 ~= 0)
    %% 4.1 execute the refine procedure for full marked(n_eg_marked = 3) with
    % style 1->4, which is illustrated as follows:
    %                               V1
    %                               /   \
    %                              /     \
    %                             /        \
    %                          (5)       (4)
    %                          /    T43     \
    %                         /                 \
    %                       V6-----(7)-----V5
    %                       /  \                  / \
    %                      /     \   T44       /   \
    %                     /        \            /      \
    %                   (6)       (8)     (9)      (3)
    %                   /   T41      \     /   T42   \
    %                  /                 \  /               \
    %               V2-----(1)------V4----(2)----V3
    tri3 = find(n_eg_marked == 3);
    % A: generate T_new information for 1->4
    T4help = [T(tri3,:)  edge_cache(TE(tri3,1),3) ...
        edge_cache(TE(tri3,2),3) ...
        edge_cache(TE(tri3,3),3)];  % cache V1 -> V6 for each 1->4 triangles
    T44 = tri3;                             Tnew(T44,:) = T4help(:,[4 5 6]);  % the center triangle replace the original one
    T41 = ((nT+1):3:(nT+3*n3))'; Tnew(T41,:) = T4help(:,[2 4 6]);
    T42 = T41 + 1                     ; Tnew(T42,:) = T4help(:,[3 5 4]);
    T43 = T41 + 2                     ; Tnew(T43,:) = T4help(:,[1 6 5]);
    
    % B: genetate TE when 1-> 4
    TE4help = zeros(n3, 9);  % unsorted subedge info
    TE4help(:,1:2) = edge_cache(TE(tri3,1), [1 2]);
    TE4help(:,3:4) = edge_cache(TE(tri3,2), [1 2]);
    TE4help(:,5:6) = edge_cache(TE(tri3,3), [1 2]);
    
    % adjust index of  subedge 1 and subedge 2
    idxNeedChange12 = ((T(tri3,2)==Enew(TE4help(:,2),1)));
    idxNeedNoChange12 = ((T(tri3,2)==Enew(TE4help(:,1),1)));
    
    tmp = TE4help(idxNeedChange12,1);
    TE4help(idxNeedChange12,1) = TE4help(idxNeedChange12, 2);
    TE4help(idxNeedChange12,2) = tmp;
    
    ETnew(TE4help(idxNeedNoChange12, 1), 1) = T41(idxNeedNoChange12);  % left triangle
    ETnew(TE4help(idxNeedNoChange12, 2), 1) = T42(idxNeedNoChange12);
    ETnew(TE4help(idxNeedChange12,     1), 2) = T41(idxNeedChange12);  % right one, eg1 and eg 2 changed
    ETnew(TE4help(idxNeedChange12,     2), 2) = T42(idxNeedChange12);
    
    % adjust index of  subedge 3 and subedge 4
    idxNeedChange34 = ((T(tri3,3)==Enew(TE4help(:,4),1)));
    idxNeedNoChange34 = ((T(tri3,3)==Enew(TE4help(:,3),1)));
    
    tmp = TE4help(idxNeedChange34,3);
    TE4help(idxNeedChange34,3) = TE4help(idxNeedChange34, 4);
    TE4help(idxNeedChange34,4) = tmp;
    
    ETnew(TE4help(idxNeedNoChange34, 3), 1) = T42(idxNeedNoChange34);
    ETnew(TE4help(idxNeedNoChange34, 4), 1) = T43(idxNeedNoChange34); %left triangle
    ETnew(TE4help(idxNeedChange34,     3), 2) = T42(idxNeedChange34);
    ETnew(TE4help(idxNeedChange34,     4), 2) = T43(idxNeedChange34);  % right triangle
    
    % adjust index of  subedge 5 and subedge 6
    idxNeedChange56 = ((T(tri3,1)==Enew(TE4help(:,6),1)));
    idxNeedNoChange56 = ((T(tri3,1)==Enew(TE4help(:,5),1)));
    
    tmp = TE4help(idxNeedChange56,5);
    TE4help(idxNeedChange56,5) = TE4help(idxNeedChange56, 6);
    TE4help(idxNeedChange56,6) = tmp;
    
    ETnew(TE4help(idxNeedNoChange56, 5), 1) = T43(idxNeedNoChange56);  % left triangle
    ETnew(TE4help(idxNeedNoChange56, 6), 1) = T41(idxNeedNoChange56);
    ETnew(TE4help(idxNeedChange56,     5), 2) = T43(idxNeedChange56);  % right triangle
    ETnew(TE4help(idxNeedChange56,     6), 2) = T41(idxNeedChange56);
    
    % new inner edge for each 1->4 triangles
    eg7 = ((nE + ne_mark + 1):3:(nE + ne_mark + 3*n3))';
    Enew(eg7,:) = T4help(:,[5 6]);   TE4help(:,7) = eg7;
    ETnew(eg7,:) = [T43 T44];
    eg8 = eg7 + 1;
    Enew(eg8,:) = T4help(:,[6 4]);   TE4help(:,8) = eg8;
    ETnew(eg8,:) = [T41 T44];
    eg9 = eg7 + 2;
    Enew(eg9,:) = T4help(:,[4 5]);   TE4help(:,9) = eg9;
    ETnew(eg9,:) = [T42 T44];
    
    % generate TE information
    TEnew(T44,:) = TE4help(:,[7 8 9]);  % the center triangle replace the original one
    TEnew(T41,:) = TE4help(:,[8 6 1]);
    TEnew(T42,:) = TE4help(:,[9 2 3]);
    TEnew(T43,:) = TE4help(:,[7 4 5]);
end

%%%%%%%
%% 4.2 execute the refine procedure for full marked(n_eg_marked = 2) with
% style 1->3, which is illustrated as follows:
%                               V1
%                               /   \
%                              /     \
%                             /        \
%                          (5)        (4)
%                          /     T31    \
%                         /                  \
%                     V5 -----(7)------ V4
%                       /                 /      \
%                      /   T32     /            \
%                   (6)           /                \
%                    /       (8)                   (3)
%                   /      /           T33          \
%                  /   /                                 \
%               V2------------(1)---------------V3
%
if ( n2 ~= 0)
        
    tri2 = find(n_eg_marked == 2);
    % !!!! classify tri1 with local index of the refined edge !!!!
    flag1 = (edge_mark(TE(tri2, 1))==0);  tri21 = tri2(flag1); % logical index is 1
    flag2 = (edge_mark(TE(tri2, 2))==0);  tri22 = tri2(flag2); % logical index is 2
    flag3 = (edge_mark(TE(tri2, 3))==0);  tri23 = tri2(flag3); % logical index is 3
    
    % A: generate T_new information
    T2help = [T(tri21,[1 2 3])  edge_cache(TE(tri21, 2), 3)  edge_cache(TE(tri21, 3), 3) ; ...
        T(tri22,[2 3 1])  edge_cache(TE(tri22, 3), 3)  edge_cache(TE(tri22, 1), 3) ; ...
        T(tri23,[3 1 2])  edge_cache(TE(tri23, 1), 3)  edge_cache(TE(tri23, 2), 3) ];  % cache V1 -> V5
    
    %tri2 = [tri21; tri22; tri23];   % reorder the tri2
    T31 = [tri21; tri22; tri23];                  Tnew(T31,:) = T2help(:, [1 5 4]);  % the center triangle replace the original one
    T32 = nT + 3*n3 + (1:2:2*n2)';  Tnew(T32,:) = T2help(:, [2 4 5]);
    T33 = T32 + 1;                                 Tnew(T33,:) = T2help(:, [3 4 2]);
    
    % B: genetate TE
    TE2help = zeros(n2, 8);  % unsorted subedge info
    TE2help(:,1) = [TE(tri21, 1); TE(tri22, 2); TE(tri23, 3) ];
    TE2help(:,[3 4 5 6]) = [ edge_cache(TE(tri21,2), [1 2]), edge_cache(TE(tri21,3), [1 2]) ; ...
        edge_cache(TE(tri22,3), [1 2]), edge_cache(TE(tri22,1),  [1 2]) ; ...
        edge_cache(TE(tri23,1), [1 2]), edge_cache(TE(tri23,2),  [1 2])  ];
%     TE2help(:,[2 7 8]) = 0;
%     
    % adjust index of  subedge 3 and subedge 4
    idxNeedChange34 = (Tnew(T33,1)==Enew(TE2help(:,4),1));
    idxNeedNoChange34 = (Tnew(T33,1)==Enew(TE2help(:,3),1));
    
    tmp = TE2help(idxNeedChange34,3);
    TE2help(idxNeedChange34,3) = TE2help(idxNeedChange34, 4);
    TE2help(idxNeedChange34,4) = tmp;
    
    ETnew(TE2help(idxNeedNoChange34, 3), 1) = T33(idxNeedNoChange34);
    ETnew(TE2help(idxNeedNoChange34, 4), 1) = T31(idxNeedNoChange34); %left triangle
    ETnew(TE2help(idxNeedChange34,     3), 2) = T33(idxNeedChange34);
    ETnew(TE2help(idxNeedChange34,     4), 2) = T31(idxNeedChange34);  % right triangle
    
    % adjust index of  subedge 5 and subedge 6
    idxNeedChange56 = (Tnew(T31,1)==Enew(TE2help(:,6),1));
    idxNeedNoChange56 = (Tnew(T31,1)==Enew(TE2help(:,5),1));
    
    tmp = TE2help(idxNeedChange56,5);
    TE2help(idxNeedChange56,5) = TE2help(idxNeedChange56, 6);
    TE2help(idxNeedChange56,6) = tmp;
    
    ETnew(TE2help(idxNeedNoChange56, 5), 1) = T31(idxNeedNoChange56);  % left triangle
    ETnew(TE2help(idxNeedNoChange56, 6), 1) = T32(idxNeedNoChange56);
    ETnew(TE2help(idxNeedChange56,     5), 2) = T31(idxNeedChange56);  % right triangle
    ETnew(TE2help(idxNeedChange56,     6), 2) = T32(idxNeedChange56);
    
    % new inner edge for each 1->4 triangles
    eg7 = (nE + ne_mark + 3*n3 + (1:2:2*n2))';
    Enew(eg7,:) = T2help(:,[4 5]);
    TE2help(:,7) = eg7;
    ETnew(eg7,:) = [T31 T32];
    eg8 = eg7 + 1;
    Enew(eg8,:) = T2help(:,[2 4]);
    TE2help(:,8) = eg8;
    ETnew(eg8,:) = [T32 T33];
    
    % UPDATE ET INFO FOR EDGE e3 in the above figure
    ldx = ((ET(TE2help(:,1),1) - T31) == 0); ETnew(TE2help(ldx,1),1) = T33(ldx);
    rdx = ((ET(TE2help(:,1),2) - T31) == 0); ETnew(TE2help(rdx,1),2) = T33(rdx);
    
    % generate TE information
    TEnew(T31,:) = TE2help(:, [7 4 5]);
    TEnew(T32,:) = TE2help(:, [7 6 8]);
    TEnew(T33,:) = TE2help(:, [8 1 3]);
end

%%%%%%%%%%%
%% 4.3  refine case 2, style 1->2 :
%
%                                V1
%                               /  |  \
%                             /    |     \
%                          /       |        \
%                      (4)        |        (3)
%                      /         (5)           \
%                    /    T21   |    T22    \
%                  /               |                \
%               V2---(1)----V4----(2)----V3
%
if (n1 ~= 0)
    tri1 = find(n_eg_marked == 1);
    % !!!! classify tri1 with local index of the refined edge !!!!
    flag1 = (edge_mark(TE(tri1,1))==1);  tri11 = tri1(flag1); % logical index is 1
    flag2 = (edge_mark(TE(tri1,2))==1);  tri12 = tri1(flag2); % logical index is 2
    flag3 = (edge_mark(TE(tri1,3))==1);  tri13 = tri1(flag3); % logical index is 3
    
    % A: generate T_new information
    T1help = [T(tri11,[1 2 3])  edge_cache(TE(tri11,1),3); ...
        T(tri12,[2 3 1])  edge_cache(TE(tri12,2),3); ...
        T(tri13,[3 1 2])  edge_cache(TE(tri13,3),3) ];  % cache V1 -> V4
    
    %tri1 = [tri11; tri12; tri13];   % reorder the tri1
    T21 = [tri11; tri12; tri13];       Tnew(T21,:) = T1help(:,[1 2 4]);  % the center triangle replace the original one
    T22 = nT + 3*n3 + 2*n2 + (1:n1)';  Tnew(T22,:) = T1help(:,[1 4 3]);
    
    % B: genetate TE
    TE1help = zeros(n1, 5);  % unsorted subedge info
    TE1help(:,[1 2]) = [ edge_cache(TE(tri11,1),[1 2]); ...
        edge_cache(TE(tri12,2),[1 2]); ...
        edge_cache(TE(tri13,3),[1 2]) ];
    TE1help(:,[3 4]) = [TE(tri11, [2 3]); ...
        TE(tri12, [3 1]); ...
        TE(tri13, [1 2]) ];
    
    % adjust index of  subedge 1 and subedge 2
    idxNeedChange12 = (Tnew(T21,2) == Enew(TE1help(:,2),1));
    idxNeedNoChange12 = (Tnew(T21,2) == Enew(TE1help(:,1),1));
    
    tmp = TE1help(idxNeedChange12,1);
    TE1help(idxNeedChange12,1) = TE1help(idxNeedChange12, 2);
    TE1help(idxNeedChange12,2) = tmp;
    
    ETnew(TE1help(idxNeedNoChange12, 1), 1) = T21(idxNeedNoChange12);
    ETnew(TE1help(idxNeedNoChange12, 2), 1) = T22(idxNeedNoChange12);
    ETnew(TE1help(idxNeedChange12,     1), 2) = T21(idxNeedChange12);
    ETnew(TE1help(idxNeedChange12,     2), 2) = T22(idxNeedChange12);
    
    % new inner edge for 1->2
    eg5 = nE + ne_mark + 3*n3 + 2*n2 + (1:n1)';
    Enew(eg5,:) = T1help(:,[1 4]);   TE1help(:,5) = eg5;
    ETnew(eg5,:) = [T21 T22];
    
    % UPDATE ET INFO FOR EDGE e3 in the above figure
    ldx = ((ET(TE1help(:,3),1) - T21) == 0); ETnew(TE1help(ldx,3),1) = T22(ldx);
    rdx = ((ET(TE1help(:,3),2) - T21) == 0); ETnew(TE1help(rdx,3),2) = T22(rdx);
    
    % set TE informations
    TEnew(T21,:) = TE1help(:,[1 5 4]);  % the left triangle replace the original one
    TEnew(T22,:) = TE1help(:,[2 3 5]);
    
    %% 5. At last, treat abnormal ETnew value, where the above strategy for
    %% ETnew fails (on the boundary).
    bnd = (ETnew(:,1) == 0);
    ETnew(bnd,1) = ETnew(bnd,2);  ETnew(bnd,2) = 0;
    
end


%% script to validate the mesh data
% trimesh(Tnew, Vnew(:,1), Vnew(:,2)); hold on;
% text(Vnew(:,1),Vnew(:,2),num2str((1:size(Vnew,1))'),'color','blue','FontSize',14); % mark vertex id
% ex = 0.5*(Vnew(Enew(:,1),1)+Vnew(Enew(:,2),1));
% ey = 0.5*(Vnew(Enew(:,1),2)+Vnew(Enew(:,2),2));
% text(ex,ey,num2str((1:size(Enew,1))'),'color','black');  % mark the edge idx
% ex = 0.333*(Vnew(Tnew(:,1),1) + Vnew(Tnew(:,2),1) + Vnew(Tnew(:,3),1));
% ey = 0.333*(Vnew(Tnew(:,1),2) + Vnew(Tnew(:,2),2) + Vnew(Tnew(:,3),2));
% text(ex,ey,num2str((1:size(Tnew,1))'),'color','red','FontSize',18);  % mark the vertex idx
% hold off;